As part of the PhD work of Bruno M. Guerreiro © 2024. If using this notebook, please cite the paper: https://doi.org/10.1021/acsbiomaterials.2c00075
Disclaimer: due to the changing nature of programming documentation, lab work developed and tacit knowledge in this notebook, please contact the author at bruno.guerreiro@fulbrightmail.org if something is not working properly. The code is not actively maintained.
# manage data and fit
import pandas as pd
import numpy as np
import plotly.graph_objs as go
import matplotlib.pyplot as plt
# first part with least squares
from scipy.optimize import curve_fit
# second part about ODR
from scipy.odr import ODR, Model, Data, RealData
# style and notebook integration of the plots
import seaborn as sns
%matplotlib inline
sns.set(style="whitegrid")
# Import data
%store -r water_Tnuc
%store -r fp025_Tnuc
%store -r fp05_Tnuc
%store -r chi_water
%store -r chi_fp025
%store -r chi_fp05
water = pd.DataFrame()
fp025 = pd.DataFrame()
fp05 = pd.DataFrame()
water['chi'] = np.log(pd.Series(chi_water))
fp025['chi'] = np.log(pd.Series(chi_fp025))
fp05['chi'] = np.log(pd.Series(chi_fp05))
water['temp'] = pd.Series(water_Tnuc) *-1
fp025['temp'] = pd.Series(fp025_Tnuc) *-1
fp05['temp'] = pd.Series(fp05_Tnuc) *-1
### Manually choosing what datapoints to fit
# Water
Twater_manFit = water.temp[120:len(water.temp)-10]
Chiwater_manFit = water.chi[120:len(water.chi)-10]
# Fucopol 0.25%
Tfp025_manFit = fp025.temp[30:len(fp025.temp)-25]
Chifp025_manFit = fp025.chi[30:len(fp025.chi)-25]
# Fucopol 0.5%
Tfp05_manFit = fp05.temp[15:len(fp05.temp)-15]
Chifp05_manFit = fp05.chi[15:len(fp05.chi)-15]
import math
def poisson_model(T, gamma, n):
### CONSTANTS ###
A = 2.9
beta = -2
Tm = 0
# division of constants
const = A/beta
# numerator
num = np.multiply(gamma, np.power(np.subtract(T, Tm), 1+n))
# denominator
denom = 1+n
# do division
tmp = const * (num/denom)
# exponentiate
return tmp
Chi_fit_water = []
for i in Twater_manFit:
Chi_fit_water.append(poisson_model(i, 0.0001, 4))
Chi_fit_fp025 = []
for i in Tfp025_manFit:
Chi_fit_fp025.append(poisson_model(i, 0.0001, 4))
Chi_fit_fp05 = []
for i in Tfp05_manFit:
Chi_fit_fp05.append(poisson_model(i, 0.0001, 4))
# Plot data
from plotly.subplots import make_subplots
fig = make_subplots(rows=1, cols=3)
# Experimental data
fig.add_trace(go.Scatter(name='Water', x=water.temp, y=water.chi, mode='markers', marker_color='rgba(0,0,0,.2)'), row=1, col=1)
fig.add_trace(go.Scatter(name='FP 0.25%', x=fp025.temp, y=fp025.chi, mode='markers', marker_color='rgba(0,0,255,.2)'), row=1, col=2)
fig.add_trace(go.Scatter(name='FP 0.5%', x=fp05.temp, y=fp05.chi, mode='markers', marker_color='rgba(255,0,0,.2)'), row=1, col=3)
# First-approximation fittings
fig.add_trace(go.Scatter(name='Water', x=Twater_manFit, y=Chi_fit_water, mode='lines', line_color='rgba(0,0,0,1)'), row=1, col=1)
fig.add_trace(go.Scatter(name='FP 0.25%', x=Tfp025_manFit, y=Chi_fit_fp025, mode='lines', line_color='rgba(0,0,255,1)'), row=1, col=2)
fig.add_trace(go.Scatter(name='FP 0.5%', x=Tfp05_manFit, y=Chi_fit_fp05, mode='lines', line_color='rgba(255,0,0,1)'), row=1, col=3)
fig.update_layout(template='simple_white', height=400, showlegend=False,
xaxis2_title='<b>-Tnuc (ºC)<b>',
yaxis_title='<b>ln(Chi)<b>',
font_family='Arial', font_color='black', font_size=16
)
fig.show()
RMSE_water = []
RMSE_fp025 = []
RMSE_fp05 = []
for i, j in zip(water.chi, Chi_fit_water):
RMSE_water.append(np.sqrt(np.mean((i-j)**2)))
for i, j in zip(fp025.chi, Chi_fit_fp025):
RMSE_fp025.append(np.sqrt(np.mean((i-j)**2)))
for i, j in zip(fp05.chi, Chi_fit_fp05):
RMSE_fp05.append(np.sqrt(np.mean((i-j)**2)))
from scipy.optimize import curve_fit
popt_water, pcov_water = curve_fit(f=poisson_model, xdata=Twater_manFit, ydata=Chiwater_manFit, p0=(0.01, 1), sigma=RMSE_water, maxfev=12000)
popt_fp025, pcov_fp025 = curve_fit(f=poisson_model, xdata=Tfp025_manFit, ydata=Chifp025_manFit, p0=(0.000001, 10), sigma=RMSE_fp025, maxfev=12000)
popt_fp05, pcov_fp05 = curve_fit(f=poisson_model, xdata=Tfp05_manFit, ydata=Chifp05_manFit, p0=(0.000001, 10), sigma=RMSE_fp05, maxfev=12000)
print('DATA FOR WATER:','\nOptimized value for gamma: ',popt_water[0], '\nOptimized value for n: ',popt_water[1])
DATA FOR WATER: Optimized value for gamma: 6.4199715065082686e-12 Optimized value for n: 10.269109653531459
print('DATA FOR 0.25% FUCOPOL:','\nOptimized value for gamma: ',popt_fp025[0], '\nOptimized value for n: ',popt_fp025[1])
DATA FOR 0.25% FUCOPOL: Optimized value for gamma: 2.7893101596601606e-15 Optimized value for n: 14.199805211490776
print('DATA FOR 0.5% FUCOPOL:','\nOptimized value for gamma: ',popt_fp05[0], '\nOptimized value for n: ',popt_fp05[1])
DATA FOR 0.5% FUCOPOL: Optimized value for gamma: 3.7948659545168766e-35 Optimized value for n: 33.37507578001123
# compute a standard deviation error from pcov
gamma_opt_water, n_opt_water = popt_water
gamma_opt_fp025, n_opt_fp025 = popt_fp025
gamma_opt_fp05, n_opt_fp05 = popt_fp05
# Water
perr_water = np.sqrt(np.diag(pcov_water))
Dgamma_water, Dn_water = perr_water
# FucoPol 0.25%
perr_fp025 = np.sqrt(np.diag(pcov_fp025))
Dgamma_fp025, Dn_fp025 = perr_fp025
# FucoPol 0.5%
perr_fp05 = np.sqrt(np.diag(pcov_fp05))
Dgamma_fp05, Dn_fp05 = perr_fp05
print("ERROR DATA FOR WATER:")
print("Gamma = %6.30f +/- %4.30f" % (gamma_opt_water, Dgamma_water))
print("n = %6.3f +/- %4.3f" % (n_opt_water, Dn_water))
print("")
print("ERROR DATA FOR FUCOPOL 0.25%:")
print("Gamma = %6.30f +/- %4.30f" % (gamma_opt_fp025, Dgamma_fp025))
print("n = %6.3f +/- %4.3f" % (n_opt_fp025, Dn_fp025))
print("")
print("ERROR DATA FOR FUCOPOL 0.5%:")
print("Gamma = %6.40f +/- %4.40f" % (gamma_opt_fp05, Dgamma_fp05))
print("n = %6.3f +/- %4.3f" % (n_opt_fp05, Dn_fp05))
ERROR DATA FOR WATER: Gamma = 0.000000000006419971506508268586 +/- 0.000000000002543765378348991293 n = 10.269 +/- 0.169 ERROR DATA FOR FUCOPOL 0.25%: Gamma = 0.000000000000002789310159660161 +/- 0.000000000000000624509598625161 n = 14.200 +/- 0.095 ERROR DATA FOR FUCOPOL 0.5%: Gamma = 0.0000000000000000000000000000000000379487 +/- 0.0000000000000000000000000000000000404543 n = 33.375 +/- 0.446
R2_water = (np.corrcoef(Twater_manFit, Chiwater_manFit)[0,1])**2
R2_fp025 = (np.corrcoef(Tfp025_manFit, Chifp025_manFit)[0,1])**2
R2_fp05 = (np.corrcoef(Tfp05_manFit, Chifp05_manFit)[0,1])**2
print('R2 of water: ', round(R2_water,3))
print('R2 of FP 0.25%: ', round(R2_fp025,3))
print('R2 of FP 0.5%: ', round(R2_fp05,3))
R2 of water: 0.784 R2 of FP 0.25%: 0.841 R2 of FP 0.5%: 0.918
# Plot data with optimized model
fig = go.Figure()
fig = make_subplots(rows=1, cols=3)
# Water
fig.add_trace(go.Scatter(name='Exp', x=water.temp, y=water.chi, mode='markers', marker_color='rgba(0,0,0,.2)'), row=1, col=1)
fig.add_trace(go.Scatter(name='Fit', x=Twater_manFit, y=poisson_model(Twater_manFit, gamma_opt_water, n_opt_water), mode='lines', line_width=15, line_color='rgba(0,0,0,.5)'), row=1, col=1)
fig.add_trace(go.Scatter(name='Fit', x=water.temp, y=poisson_model(water.temp, gamma_opt_water, n_opt_water), mode='lines', line_color='rgba(0,0,0,.8)'), row=1, col=1)
# FP 0.25%
fig.add_trace(go.Scatter(name='Exp', x=fp025.temp, y=fp025.chi, mode='markers', marker_color='rgba(0,0,255,.2)'), row=1, col=2)
fig.add_trace(go.Scatter(name='Fit', x=Tfp025_manFit, y=poisson_model(Tfp025_manFit, gamma_opt_fp025, n_opt_fp025), mode='lines', line_width=15, line_color='rgba(0,0,255,.5)'), row=1, col=2)
fig.add_trace(go.Scatter(name='Fit', x=fp025.temp, y=poisson_model(fp025.temp, gamma_opt_fp025, n_opt_fp025), mode='lines', line_color='rgba(0,0,255,.8)'), row=1, col=2)
# FP 0.5%
fig.add_trace(go.Scatter(name='Exp', x=fp05.temp, y=fp05.chi, mode='markers', marker_color='rgba(255,0,0,.2)'), row=1, col=3)
fig.add_trace(go.Scatter(name='Fit', x=Tfp05_manFit, y=poisson_model(Tfp05_manFit, gamma_opt_fp05, n_opt_fp05), mode='lines', line_width=15, line_color='rgba(255,0,0,.5)'), row=1, col=3)
fig.add_trace(go.Scatter(name='Fit', x=fp05.temp, y=poisson_model(fp05.temp, gamma_opt_fp05, n_opt_fp05), mode='lines', line_color='rgba(255,0,0,.8)'), row=1, col=3)
fig.update_layout(template='simple_white', showlegend=False, height=400,
xaxis1_title='',
xaxis3_title='',
xaxis2_title='<b>-Tnuc (ºC)<b>',
xaxis1_range=[4,15],
yaxis_title='<b>ln(Chi)<b>',
font_family='Arial', font_color='black', font_size=16
)
fig.show()
fig = make_subplots(rows=1, cols=3)
# Water
fig.add_trace(go.Scatter(name='Exp', x=np.exp(water.chi), y=water.temp*-1, mode='markers',marker_color='rgba(0,0,0,.2)', marker_size=8, marker_symbol=100),row=1,col=1)
fig.add_trace(go.Scatter(name='Fit', x=np.exp(poisson_model(water.temp, gamma_opt_water, n_opt_water)), y=water.temp*-1, mode='lines', line_color='rgba(0,0,0,1)'),row=1,col=1)
# FP 0.25%
fig.add_trace(go.Scatter(name='Exp', x=np.exp(fp025.chi), y=fp025.temp*-1, mode='markers',marker_color='rgba(0,0,255,.2)', marker_size=8, marker_symbol=100),row=1,col=2)
fig.add_trace(go.Scatter(name='Fit', x=np.exp(poisson_model(Tfp025_manFit, gamma_opt_fp025, n_opt_fp025)), y=Tfp025_manFit*-1, mode='lines', line_color='rgba(0,0,255,1)'),row=1,col=2)
# FP 0.5%
fig.add_trace(go.Scatter(name='Exp', x=np.exp(fp05.chi), y=fp05.temp*-1, mode='markers',marker_color='rgba(255,0,0,.2)', marker_size=8, marker_symbol=100),row=1,col=3)
fig.add_trace(go.Scatter(name='Fit', x=np.exp(poisson_model(fp05.temp, gamma_opt_fp05, n_opt_fp05)), y=fp05.temp*-1, mode='lines', line_color='rgba(255,0,0,1)'),row=1,col=3)
fig.update_layout(template='simple_white', width=1000, height=500, showlegend=False,
xaxis2_title='<b>Unfrozen fraction (%)<b>',
yaxis_title='<b>Nucleation temperature (ºC)<b>',
yaxis1_range=[-15,-5],
font_family='Arial', font_color='black', font_size=16,
xaxis_range=[-0.05,1.05])
fig.update_xaxes(mirror=True, showline=True, ticks='outside', linewidth=2)
fig.update_xaxes(mirror=True, showline=True, ticks='outside', linewidth=2, tickprefix="<b>",ticksuffix ="</b>")
fig.update_yaxes(mirror=True, showline=True, ticks='outside', linewidth=2, dtick=1, tickprefix="<b>",ticksuffix ="</b>")
fig.show()
fig = go.Figure()
# Water
fig.add_trace(go.Scatter(name='Exp', x=np.exp(water.chi), y=water.temp*-1, mode='markers',marker_color='rgba(0,0,0,.8)', marker_size=8, marker_symbol=115))
#fig.add_trace(go.Scatter(name='Fit', x=np.exp(poisson_model(water.temp, gamma_opt_water, n_opt_water)), y=water.temp*-1, mode='lines', line_color='rgba(0,0,0,.25)'))
# FP 0.25%
fig.add_trace(go.Scatter(name='Exp', x=np.exp(fp025.chi), y=fp025.temp*-1, mode='markers',marker_color='rgba(0,0,255,.8)', marker_size=8, marker_symbol=15))
#fig.add_trace(go.Scatter(name='Fit', x=np.exp(poisson_model(fp025.temp, gamma_opt_fp025, n_opt_fp025)), y=fp025.temp*-1, mode='lines', line_color='rgba(0,0,255,.25)'))
# FP 0.5%
fig.add_trace(go.Scatter(name='Exp', x=np.exp(fp05.chi), y=fp05.temp*-1, mode='markers',marker_color='rgba(255,0,0,.8)', marker_size=8, marker_symbol=15))
#fig.add_trace(go.Scatter(name='Fit', x=np.exp(poisson_model(fp05.temp, gamma_opt_fp05, n_opt_fp05)), y=fp05.temp*-1, mode='lines', line_color='rgba(255,0,0,.25)'))
fig.update_layout(template='simple_white', width=350, height=500, showlegend=False,
xaxis_title='<b>Unfrozen fraction (%)<b>',
yaxis_title='<b>Nucleation temperature (ºC)<b>',
font_family='Arial', font_color='black', font_size=14,
xaxis_range=[-0.05,1.05])
fig.update_traces(marker_size=8, marker_symbol=0)
fig.update_xaxes(mirror=True, showline=True, ticks='outside', linewidth=2)
fig.update_yaxes(mirror=True, showline=True, ticks='outside', linewidth=2, dtick=1)
fig.update_xaxes(mirror=True, showline=True, ticks='outside', linewidth=2, tickprefix="<b>",ticksuffix ="</b>")
fig.update_yaxes(mirror=True, showline=True, ticks='outside', linewidth=2, dtick=1, tickprefix="<b>",ticksuffix ="</b>")
fig.show()
#fig.write_image("images/survivor.svg")
J_water = np.multiply(gamma_opt_water, np.power(water.temp, n_opt_water))
J_fp025 = np.multiply(gamma_opt_fp025, np.power(fp025.temp, n_opt_fp025))
J_fp05 = np.multiply(gamma_opt_fp05, np.power(fp05.temp, n_opt_fp05))
# Fitting
T_range = np.linspace(0.01,16, 100)
IT_water_fit = 1/np.multiply(gamma_opt_water, np.power(T_range, n_opt_water))
IT_fp025_fit = 1/np.multiply(gamma_opt_fp025, np.power(T_range, n_opt_fp025))
IT_fp05_fit = 1/np.multiply(gamma_opt_fp05, np.power(T_range, n_opt_fp05))
IT_water = 1/J_water
IT_fp025 = 1/J_fp025
IT_fp05 = 1/J_fp05
#####
#IT_waterOptimal_fit = 1/np.multiply(gamma_opt_water, np.power(T_range, n_opt_water))
#IT_waterRepeat_fit = 1/np.multiply(6.23e-12, np.power(T_range, 7.8))
fig = go.Figure()
# Fits
fig.add_trace(go.Scatter(name='Water', x=T_range*-1, y=IT_water_fit, mode='lines', line_width=3, line_color='rgba(0,0,0,.75)'))
fig.add_trace(go.Scatter(name='FucoPol 0.5%', x=T_range*-1, y=IT_fp05_fit, mode='lines', line_width=3, line_color='rgba(255,0,0,.75)'))
fig.add_trace(go.Scatter(name='FucoPol 0.25%', x=T_range*-1, y=IT_fp025_fit, mode='lines', line_width=3, line_color='rgba(0,0,255,.75)'))
fig.update_layout(template='simple_white', width=600, height=500, showlegend=True,
yaxis_title='<b>Induction Time<b>',
xaxis_title='<b>Working temperature (ºC)<b>',
font_family='Arial', font_color='black', font_size=14)
tickvals = [1/60000, 1/60, 1, 60, 1440, 10080, 40320, 483840, 4838400, 48384000, 483840000, 48384000000]
ticktext = ['Millisecond','Second', 'Minute', 'Hour', 'Day', 'Week', 'Month', 'Year', 'Decade', 'Century', 'Millennium', '∞']
fig.update_xaxes(mirror=True, showline=True, ticks='outside', linewidth=2, tickprefix="<b>",ticksuffix ="</b>",
range=[-4, -13])
fig.update_yaxes(mirror=True, showline=True, ticks='outside', linewidth=2, dtick=1, type='log',
tickprefix="<b>",ticksuffix ="</b>", tickvals=tickvals, ticktext=ticktext, range=[-0.5,7])
fig.show()